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We present an exact treatment of the hysteresis behavior of the zero-temperature random-field 
Ising model on a Bethe lattice when it is driven by an external field and evolved according to a 
2-spin-flip dynamics. We focus on lattice connectivities z = 2 (the one-dimensional chain) and 
2 = 3. For the latter case, we demonstrate the existence of an out-of-equilibrium phase transition, 
in contrast with the situation found with the standard 1-spin-flip dynamics. We discuss the influence 
of the degree of cooperativity of the (local) spin dynamics of the nonequilibrium response on the 
system. 
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I. INTRODUCTION 

The nonequilibrium response of the random field Ising model (RFIM) at zero temperature to a smoothly chang- 
ing external magnetic field has been introduced by Sethna et al.P, 01 to describe the hysteresis behavior of low- 
temperature magnets, in particular the property of "return-point memory" and the power-law scaling characteristic 
of the "Barkhausen noise" 0. In such a description, the competition between the quenched disorder and the fer- 
romagnetic interactions induces a very large number of local energy minima in which the system gets trapped. In 
the zero-temperature dynamics, no escape by thermal activation is possible and evolution only results from changes 
in the driving field. Instead of wandering to find energetically more favorable states, the system then goes to the 
closest minimum made accessible by the field change. The state of the system is history-dependent and the evo- 
lution proceeds by "avalanches", i. e., jumps from one mini! mum to another. An extensive investigation of this 
zero-temperature hysteretic behavior of the RFIM has been carried out by several authors (for a recent review, see 
0). It has been shown that the model displays an out-of-equilibrium critical point, with its associated scaling and 
universality properties. The critical point separates a strong-disorder regime in which the magnetization curve and 
the corresponding hysteresis loop is smooth on the macroscopic scale and a weak-disorder regime in which it shows a 
discontinuous jump at a nonzero value of the external field. 

Real materials, however, are not at zero temperature. Application of the above picture requires a separation of 
time scales: the time for "local equilibration", namely the duration of the avalanche that follows a change in the 
applied field, must be much shorter than the experimental time scale, set by the rate of change of the applied field, 
which itself must be much less than the characteristic lifetime of the relevant metastable states in which the system is 
trapped. The extreme limit of the former condition is the "adiabatic" case in which the rate of change of the external 
field goes to zero. The extreme limit of the latter condition corresponds to the zero-temperature (T = 0) dynamics in 
which no thermally activated escape takes place. The effect of a finite rate has been fairly extensively studied[^|^; 
that of a finite temperature is a tougher issue. 

To study the robustness of the scenario found for the T — Q RFIM with the standard single-spin-flip dynamics, it 
has been proposed to consider the nonequilibrium response of the RFIM still at T = 0, but with dynamics involving 
possible cooperative flips of several spinsf^. In such dynamics, the system can effectively overcome energy barriers 
that would block it in a one-spin- flip dynamics. Even if thermally activated processes are absent, the metastable states 
visited by the system when the applied field is ramped up or down are a priori different from those reached in the 
one-spin-flip dynamics. A simulation of the T = RFIM with 2-spin-flip dynamics on a cubic lattice showed that the 
change of dynamics yields a significant reduction of the width of the hysteresis loop (or, in magnetic language, of the 
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coercivity), but that the out-of-equihbrium critical behavior is still present and appears, within statistical uncertainty, 
to belong to the same universality c! lass as that observed when using the standard one-spin-flip dynamics |7)|. 

In this article, we study the hysteresis behavior of the RFIM on a Bethe lattice with a 2-spin-flip dynamics at zero 
temperature. We consider a gaussian distribution of the random fields. Our goal is to obtain an exact solution for 
the main hysteresis loop as a function of disorder strength, and to compare the results with those of the 1-spin-flip 
dynamics already derived^, 0|- A question that was found intriguing in the latter case is the absence of out-of- 
equilibrium transition for a Bethe lattice with connectivity z=3[9|. A transition is found for z > 4 and no transition 
occurs when z = 2 (the one-dimensional case). The fact that a qualitative change takes place for z > 3 and not, 
as found in equilibrium, for z > 2 is, to the least, unusual. It casts a doubt on the conjecture that the out-of- 
equilibrium and equilibrium critical points belong to the same universality class and are controlled by the same fixed 
point |1. Hol ITU . ! We thus focus on the Bethe lattice with connectivity z = 3 and we show that an out-of-equilibrium 
transition does occur in this case with the 2-spin-flip dynamics. 

The rest of the paper is organized as follows. In section II we present the model and the dynamics, contrasting the 
2-spin-flip dynamics with the standard 1-spin-flip one and reviewing the main properties of the 2-spin-flip dynamics. 
The next section is devoted to deriving the self-consistent equations that allow to calculate the main hysteresis loop 
on a Bethe lattice. In section IV, we consider the one-dimensional chain (z = 2) and in section V, the z = 3 case. 
In both cases, we compare with the results obtained with the 1-spin-flip dynamics as well as with the equilibrium 
(ground state) behavior. For the connectivity z = 3, we show that, in contrast with the 1-spin-flip dynamics, the 
2-spin-flip dynamics leads to a disorder-induced, out-of-equilibrium transition. Finally, in the last sections, we discuss 
possible extensions of the approach to the calculation of other observables and to more general k-spin-flip dynamics, 
with an emphasis on the influence of the cooperativ! ity of the (local) dynamics on the hysteresis properties, and we 
conclude. 



II. MODEL AND DYNAMICS 



We study the RFIM on a lattice of connectivity z, with Ising spin variables Si — ±1 ferromagnetically interacting 
with their neighbors and subject to both an intrinsic random field and an applied uniform field H. The Hamiltonian 
is given by 

H=-jJ2 S^SJ - 5](/i. + H)S,, (1) 

<ij> i 

where the sum < ij > is over all distinct "bonds" of the lattice and the random fields hi on each site are independently 
drawn from a gaussian probability distribution of zero mean and variance equal to cr, p{h) = exp(— /i^/2(T^)/\/27rCT . 

We consider the limit of zero temperature and the evolution of the system as it is (adiabatically) driven by the 
external magnetic field H . At T = 0, metastability results from the absence of full equilibration and the presence of 
several distinct configurations in which the system can be trapped. It is customary to define such configurations as the 
local energy minima, namely the states that minimize the Hamiltonian in Eq. (1) (due to the presence of the applied 
magnetic field, the energy is actually a magnetic enthalpy, but in the following we shall keep a loose usage of the term 
"energy"), and to envisage the phenomenon within the topographic picture of a rugged energy landscape. For Ising 
spins, the energy levels are discretized and the topographic view requires some adjustments. A useful characterization 
of the minima consists in sorting them according to the number of spins that must be flipped to yield a lower-energy 
state. More precisely, one can define a fc-stab! le minimum as a configuration whose energy cannot be decreased by 
flipping any subset of fc, or less, spins [T^. The usual minima are then 1-stable configurations, which are stable with 
respect to the flip of any single spin in the system. The ground state of the system on the other hand corresponds 
to a configuration that is stable with respect to any arbitrary large number of (simultaneous) spin flips, i. e., to the 
limit k ooOIll. 

For classical spins with no intrinsic dynamics, knowledge of the energy landscape is not enough if one does not 
provide a dynamical rule to explore the landscape. For the problem of interest here, most studies have considered the 
standard 1-spin-flip dynamics. Under such dynamics, the system evolves along irreversible paths among local energy 
minima which by construction are 1-stable for the relevant range of applied field H. At each site i, the spin Si aligns 
with the net local field fi, with 

f,^ jY,Sj+h,+H, (2) 

m 



where the sum is over all connected neighbors of site i. 
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A natural extension of this dynamics consists in allowing cooperative flips of pairs of spins. The system then jumps 
from one state to another (via an avalanche) whenever the former state becomes unstable to the flip of any single 
spin or to the simultaneous flip of any single pair of spins. Evolution now proceeds among 2-stable local minima. 
New features with respect to the 1-spin-flip dynamics are only introduced when the flipping pair is made up by two 
neighboring spins, and the local energy to be minimized is then 

Hij — —fj^Si — fjSj — JSiSj (3) 

where 

/; Sk + h, + H (4) 

is the local fleld experienced by Si without the influence of the neighbor Sj (the sum in the above expression is over all 
neighbors of i except j). The transition of the pair (ij) from one state to another under the 2-spin-flip dynamical rule 
is graphically illustrated in Figure 1. One can think the dynamics as made of both single-spin-flips and "irreducibly 
cooperative" 2-spin flips. The former occur whenever the net local fleld on Si, ft = // + J Si, or Sj, fj = /j + J Si, 

changes sign. This corresponds to the transitions < — > — h , < — > H — and — h < — > ++ , H — < — > ++ in the 

diagram. An irreducibly cooperative 2-spin-flip involves a pair of neighboring spins which arc of the same sign ( 

or H — h) and cannot flip individually, i.e., without the simultaneous flip of the neighbor. This takes place whenever the 
net fleld on the pair of ahgned spins, fij = fi+lf'j, changes sign in the central region of the diagram: —J < f- < J 
(condition for forbidding spin Si to flip individually) and —J<fj<J (condition for forbidding spin Sj to flip 
individually). The applied magnetic fleld is varied until a flrst pair of spins become marginally stable. This pair can 
be represented by a point on the diagram in Figure 1 and the further evolution of the pair is dictated by the border 
that is flrst attained (notice that the axes // and /j depend linearly on H so that the displacement of the point is 
always along the (1, 1) diagonal). 



(-+) 




+J 


(++) 






+J 








// 


(— ) 


-./ 




(+-) 



FIG. 1: Stability diagram for a pair of nearest-neighbor spins i,j under the 2-spin-flip dynamical rule. The axes are the local 
fields defined in Eq. (4) . 



In Ref.j3], it was shown that, due to the ferromagnetic nature of the interactions, the 2-spin-flip dynamics shares 
with the 1-spin-flip one important properties, known as the no-passing rule, the abelian property, and the existence 
of return-point memory. As a result, (i) the flnal 2-stable conflguration reached through an avalanche is independent 
of the order in which the unstable spins have been "relaxed", i. e. flipped to minimize the energy), (ii) the state 
reached at a given applied fleld is independent of the details of the previous history and only depends on the state in 
which the fleld was last reversed, provided the fleld is adiabatically varied, and (iii) in the relaxation process between 
two metastable states after a change of the applied fleld, no spins flip more than once. 

In the following we focus on the main hysteresis loop obtained by starting from a satured conflguration with all 
spins pointing in the same direction in a very large (positive or negative) applied fleld. Because of the symmetry of 
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the model, the upper branch of the loop followed by the system when decreasing the applied field from +00 to —ao 
is easily deduced from the lower branch obtained by increasing the field from —00 to +00, and we therefore only 
consider the latter. We search for an exact solution of the problem on a Bethe lattice of connectivity z. 



III. SELF-CONSISTENT EQUATIONS ON A BETHE LATTICE 



As is standard in statistical physics, we consider a Bethe lattice as the "deep interior" of a Cayley tree in the limit 
where the number of sites in both go to infinity, the size of the whole tree going to infinity faster than that of the 
chosen interior piece p^. The tree structure allows to write recursion relations and the above trick allows to close the 
relations by removing the effect of the boundaries. A detailed analysis of the hysteresis behavior of the T = RFIM 
on a Bethe lattice with the 1-spin-flip dynamics has already been Drovided|M M. fill UA. Hfll lift Il7| . Extension to more 
general dynamics is however nontrivial and we detail here the main aspects of the formalism. 

We first choose at random a site in the deep interior of the Cayley tree of connectivity z and take it as the origin. 
From this origin, one can label the "descendant" sites according to their level in the tree. A spin at level n will then 
have one parent spin at level {n — 1) and z — 1 descendant spins at level (n + 1): see Figure 2. Our goal is to compute 
the average magnetization as a function of the magnetic field H in the lower branch of the hysteresis loop. By using 
the properties of the dynamics, it is then sufficient to start from a configuration where all spins are down and to follow 
the relaxation of the spins starting from the outer part of the lattice and progressing inward, toward the origin, until 
a 2-stable configuration at the given value H of the field is reached. 




FIG. 2: llustration of the branching characteristics of the Cayley tree for 2 = 4. 



In the 1-spin-flip dynamics, the central quantity is the conditional probability that a spin at level n is flipped in 
the relaxation process, given that its parent at level {n— 1) is kept down. This conditional probability is a function of 
the external field only (it is already an average with respect to the random fields) and it obeys a recursion equation 
whose fixed point describes the typical behavior on the Bethe lattice. Deriving the counterpart of this quantity for 
the 2-spin-flip dynamics is a little tricky. Indeed, cooperative flips of pairs of spins couple the levels in a more intricate 
way. The appropriate quantity is actually the conditional probability that a spin at level n is in state S (with = ±1, 
where as usual " and "up/down" are used interchangeably) given that: (i) its local random field is h, (ii) all 

spins (and pairs of spins) in the branch are relaxed with the constraint that the parent spin at level (n! — 1) is kept 
down, and (iii) the z — 1 descendant spins at the level {n + 1) are in states 81,82, ■■,8z-i- The probability is now 
a function of the sum of the random and the uniform fields, x = h + H, and it is easy to show that it depends on 
the state of the descendant spins at level (n + 1) only through the sum 8 = '^i- We denote this conditional 

probability p^^^ (x) . 

To write the recursion relation for this probability, one needs to introduce 2-spin-fiip transition rates involving the 
change of the configuration formed by a spin at level n and its 2; — 1 direct descendants, given the state of the parent 
at level (n — 1) and the state of all the (z — 1)^ descendants at level (n + 2). Let us denote 

^{S(„+2)|(s;„^^),s;„))^(S(„+i),S(„))|S(„_i)}(^("+i)' 2:(n)) (5) 

the transition rate for going from the configuration {^[n+i)' ^[n))^ where S'^^,^-^^ = {5( , S'2_i}(„+i) to the config- 
uration (S(„+i), ^(n)), where S(„+i) = {81, 8z-i}i^n+i), given that the parent at level {n — 1) is in state S'(„_i) 
and the descendants at level (n -|- 2) are in states S(„+2) = {Si, 8z-i}(n+2) (again, it is easy to show that the rate 
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only depends on the sum of the states of the z — 1 descendants merging at the common root at level (n + 1), e.g., 
Si = J2i~i 'S'j^i); the rate depends on the local fields at level n, a;(„), and at level (n + 1), X(„_|_i) = {xi, Xz-i}n+i- 
We shall only be interested in the situation where 5^,^^ = —1 and 6'(„_i) = —1. A graphical representation is given in 
Figure 3. The transition rate can be written by combining stabilit! y diagrams similar to that in Figure 1 for all the 
pairs of spins involved in the change of state. An explicit expression will be provided below for the case z = 2. 





FIG. 3: Schematic illustration of the transition rate W^g^ +2)l(s' s' )^(S( +i) S( ))I'S( involved in the 

recursion equations for the conditional probabilities (see Eq. (5)). The configuration formed by a spin (initially down) and its 
z — 1 descendants changes via both single-spin flips and irreducibly cooperative 2-spin flips, given that the parent (the rightmost 
spin) is kept down and that the descendants on the left are in states {Si, . . . Sz-i)- 



The recursion relation for the conditional probability (x) can now be derived by considering all ways to generate 
the chosen configuration at level n from configurations on the z — 1 branches at level (n + 1). One obtains 

x^{s<„,.,,(s;„,,,-i)^(s<„,,„5|-i}('^("+i)'^) (6) 

where we have used the short-cut notation 6x = p{x — H)dx and {Si . . . S^-i} is any configuration such that 5*1 -|- 
. . . + Sz-i = S; as already stressed, it is indeed easy to prove that the conditional probability is the same for all such 
configurations: p^"!(.xO = g sums on the 5j's are restricted to the values 2r — (z — 1) with r an 

integer between and z — 1 (the combinatorial factors then count the number of different configurations of 2; — 1 spins 
with a sum equal to 2r — (z — 1)). The recursion equation can be closed by making the "homogeneous" assumption 
that all sites in the deep interior are equivalent, provided they have the same random field. This amounts to consider 
the fixed point of the recursion, which is solution of the self-consistent equation obtained by dropping the superscripts 
n and n + 1 in the above equation. 

One more step is still required to compute the average magnetization. This latter is simply related to the average 
probability that the central spin (at the randomly chosen origin of the Bethe lattice) is up: 

m{H)=2P{+l;H)-l, (7) 
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where the overbar denotes an average over the random field distribution, P(+1:H) = J SxP{+l: x, H). The prob- 
abihty that the central spin is up is obtained by studying the merging of the z branches of the tree at the origin. 
The state of each branch given that the central spin is kept down is independent from that of the other branches 
and is characterized by the conditional probability Pgg{x). Full relaxation implies to relax the constraint on the 
central spin and let the system evolve according to the 2-spin-flip dynamical rule. This evolution is described with 
the help of transition rates W that are given by an extension of those introduced above: the transition now involves 
configurations formed by the central spin and its z neighbors, given the states of the direct descendants of those 
neighbors, see Figure 4. One then obtains 

P{+1;H)= 6x Y ^xi---Sx,pg^g^{xi)---pg^s^{x^) 

^^{Si...S;,l(s;...s^-i)^(Si...s,,+i)}(^i ■ ■ ■^z-.x), 

where W is the transition rate defined just above; to simplify the notations we have included, here and in the rest of 

the paper, the combinatorial factor cf_X'^~'^^^'^ in a redefinition of the conditional probability Pssi^)- 

Having given the general formalism to calculate the main hysteresis loop on a Bethe lattice under a 2-spin-flip 
dynamics, we now specialize the study to the cases z = 2 and z = 3. 





FIG. 4: Schematic illustration of the transition rate W^^Sj...s^|(sj...s' ,-i)^(Si...s^,+i)}i^'^ . . .XzjX) involved in the expression of 
the magnetization (see Eq. (8)). The central spin flips up in a change of configuration that involves its z neighbors, given the 
states of the descendants {Si, . . . Sz)- The change of configuration involves both single-spin fiips and irreducibly cooperative 
2-spin flips. 



IV. EXACT SOLUTION FOR THE ONE-DIMENSIONAL CHAIN 



The case of the one-dimensional chain (z = 2) is of course the easiest to handle. The self-consistent equation 
obtained from Eq. (6) is linear in the conditional probability, but in contrast with what occurs for the 1-spin-flip 
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dynamics, the probability is a nontrivial function of the local field and Eq. (6) has the form of an integral equation: 

We first specify the explicit form of the transition rates that appear in this equation. We consider the transition 

for a pair of neighbors given that the descendant is in state S and that the parent is down. Recalling the relation 
between the local fields f[ introduced in section II and the local field Xi used in section III (namely, is equal to 
Xi plus the field due to the interaction with the neighbor that does not belong to the pair under consideration) and 
that the external field H increases, we can adapt the stability diagram of Figure 1 to describe the different terms and 
sectors composing the transition rates. The resulting expressions are given below: 

W^S,\{—)^{++)\-}i.Xi,x) ^e{0<x< 2J)0{~J{Si ~l)-x<xi< -J{Si - 1))+ 

9{x > 2J)9{-JiSi + 1) < xi < -J{Si - 1)) 





(— 


)^(-+) 




x) 


= 9 


[x > 2J)9{xi < 


-J{Si 


+ 1)) 




(— 


)^(+-) 




x) 


= 










(— 


)^(— ) 




x) 


= 9 


[x < 0)9{xi < - 


-J{Si - 


l)) + 9{0 




(+- 


)^(++) 




x) 


= 9 


[x > 0)9{xi > - 


-J{Si - 


1)) 




(+- 


)^(-+) 




x) 


= 










(+- 


)-(++) 




x) 


= 9 


[x < 0)9{xi > - 


-J{Si - 


1)) 




(+- 






x) 


= 









(10) 

where 9{...) is the characteristic function of the domain of local field indicated by the argument: it is 1 in the domain 
and outside; with this definition, 6{x > 0) is the usual Heaviside function H{x), 9{0 < x < 2 J) is equal to 
H{x) — H{x — 2J). and so on. 

By inserting the above expressions into Eq. (6), one finds that the conditional probabilities satisfy the following 
equations: 



p++{x)=e{Q<x<2J) 
+ e{x > 2J) 

p-+{x) = 9{x > 2J) 
p+-{x) = 9{x < 0) 

p—{x) = e{x < 0) 



E 

. Si 



+ 0O 



J(Si-l) 



5xiPs^^{xi) + 



Si 



-./(Si-l) 



J(Si-l)-x 



SxiPs^_{xi) 



E 

Si 

E 

. Si 



+ 00 



J(Si-l) 



Sxipg^^ixi) 



Si 



-J(Si-l) 



J(Si+l) 



SxiPs^_{xi) 



E 



+ 00 



J(Si-l) 



Sxip§^_{xi) 



SXlPg^^{Xl) 



E 

. Si 



-J(Si-l) 



Sxip§^_{xi) 



+ e{0<x< 2 J) 



E 

. Si 



-J{Si-l)-x 



Sxip§^_{xi) 



(11) 



Notice that the dependence on x comes through both the characteristic functions and, in the case of p++ and p- 
through an additional dependence in the region 0<x<2J (see the bounds of the integrals). 
It is convenient to introduce the disorder-averaged values of the conditional probabilities. 



Pss = / ^^Pssi^)- 



(12) 
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FIG. 5: (Color on line) Conditional probabilities p§g{x) as a function of the local field x for the 1-dimensional chain {z — 2) 
and for a = 1.0. Different values of the applied field are considered (here and in all figures the fields and the disorder strength 
are in units J = 1). Note the continuity properties discussed in the text. 

These quantities can be computed directly in computer simulations, as will be discussed below. As detailed in 
Appendix A, the solution to the set of coupled linear integral equations formed by Eqs. (11) can be found by 
reexpressing all the terms needed in the expression of the conditional probabilities with the help of the disorder- 
averaged values and of a function q(x) that is nonzero only in the domain < a; < 2 J (where it is equal to p (x), 

by definition). After some algebra, one finds 

/■2J 

P++ =PoPi(l +P2 -Pi) +Pi -Po - [1 +Po(j32 -Pi)] / 5xq{x) 

Jo 

.2,7 

P+- =Pi(l -Pi) - (1 -Pi) / Sxq{x) 

Jo 

P-+ =Po[l -Pi(l +P2 -Pi)] +Po(P2 -Pi) / Sxq{x) 

Jo 

.2.7 

p__ = (l-pi)2 + (2-pi) / Sxqix) (13) 

Jo 

where the function q(x) satisfies an inhomogeneous integral equation, 

p2J-x pO pO p2J 

q{x)- 6xiq{xi) ^ 1 - pi{l + 6xi) + ( SxAl 5xiq{xi) (14) 



— X J —X 



for Q < X < 2J. In the above expressions, we have introduced the standard probabilities that play a key role for the 
1-spin-flip dynamics, Pn{H) — J(^°°2n)j = I{z'^2n)j-H '^^P^^^ ' with n = 0, 1, or 2: Pn{H) is the probability that a 
spin with n neighbors up and z — n neighbors down flips up in an applied field H via the 1-spin-flip mechanism. The 
probabilities Pgg{x) satisfy the normalization property 



P++{x) + p+-{x) + p-+{x) + p — [x) = 1, 



(15) 
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for any x. (As a result the normalization condition is also satisfied by the disorder-averaged probabilities, as can be 
verified by using Eqs. (13).) It is also easy to check that the functions P§g{x) satisfy continuity relations: for instance, 

p++(0+) = (0^) and p__(2J~) — p ^(2J+). Finally, one can recover the results of the 1-spin-flip dynamics0 

by setting q{x) = (1 — pi)/{l — Pi + Po) over the whole interval [0, 2J]. 




-1 -0.5 0.5 1 1.5 2 
H 



FIG. 6: (Color on line) Average conditional probabilities as a function of the applied field H for z — 2 and a — 1.0. The lines 
correspond to the exact solution of Eqs. (13,14) and the symbols to the simulation results discussed in section VI [Nsites ~ 10*, 

^realizations — 10 ). 



We have numerically obtained the solution of Eq. (14) (the solution can be computed with any arbitrary precision), 
from which we have calculated the conditional probabilities and their average. The resuts are illustrated in Figures 5 
and 6. The continuity relations are clearly seen on Figure 5. One moreover finds that the average probabilities are 
smooth (continuously differentiable) functions of the external field for all values of the disorder strength cr. 

Going from the conditional probabilities to the average magnetization needs the additional step provided by Eq. 
(8). The expressions of the transition rates are now a little more complicated than those in Eqs. (10), since one must 
consider the different transition possibilities for a configuration of 3 spins (see Figure 4 with z = 2). The explicit 
formulae are given in Appendix A 

The resulting magnetization curve m{H) corresponding to the lower branch of the hysteresis loop is shown in Figure 
7 for several values of the disorder strength a. Again, the curves are always smooth, irrespective of the value of a. As 
expected, no out-of-equilibrium transition takes place in the 2-spin-flip dynamics of the T — RFIM in one dimension. 
This will be discussed in more detail below. 




FIG. 7: (Color on line) Lower branch of the hysteresis loop for several values of the disorder {z — 2). The lines correspond to 
the exact solution and the symbols to the simulation results {Nsues = 10**, Nreaiizations ~ 10^). 
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V. BETHE LATTICE WITH CONNECTIVITY z = 3 



As mentioned in the Introduction, the case ^ = 3 is especially interesting. Whereas a disorder-induced transition 

(critical point) is found in the ground state, no such transition occurs in the noncquilibrium response with 1-spin-flip 
dynamics. The calculation of the hysteresis loop under 2-spin-flip dynamics is, however, considerably more involved 
than in the one- dimensional case. 




FIG. 8: Conditional probabilities p§g{x) as a function of the local field x for z = 3 {H = 0.4 and a = 1.1). The sum of the 
probabilities is equal to 1 in all intervals. The plateaus correspond to the constants pgg which, as the functions in the intervals 
[—J, J] and [J, 3J], are defined in Eqs. (16). 



Specializing the self-consistent equations for the conditional probabilities, Eqs. (6), to the case z = 3 leads to 
nonlinear coupled integral equations. The explicit form of the transition rates must now be derived by considering 
the two pairs formed by a spin and its two descendants and by combining two diagrams similar to that of Figure 1 to 
form a three-dimensional diagram with three local fields as axes. The number of different conditional probabilities, 
p§ si^)^ is six (corresponding to 5 = — 2,0, -|-2 and S = +, — ), because the dependence on the states of the two 
descendant spins comes only through their sum. The dependence of these conditional probabilities on the local field 
X is piecewise, as in the z = 2 case. However, contrary to the z = 2 case, several nontrivial functions of x in the 
intervals [—J, J] and [J, 3 J] must be introduced instead of only one . As discussed in more detail in Appendix B, the 
general form of! the conditional probabilities is as follows: 



P+2,+ i^' 

P-2,+ {x 
P+2-{x 

Po,-{x 

P-2-{x 



= 9{-J <x < J)r2{x) + e{J <x < SJ)q2{x) + e{x > 3J)p+2,+ 
= 0{J <x < 3J)qi{x) + e{x > 3J)po,+ 
= e{x > 3 J)p_2,+ 

= e{x < -j)p+2- 

= e{x < -J)po- + 6{-J <x< J)ri{x) 
= 6{x < J)p-2- + 0{-J <x< 3J)qo{x) 



(16) 



where the expressions for the 6 constants, Pa,a, and the 5 functions, qo{x), qi{x), q2{x) and rx{x) and r2{x), are given in 
Appendix B. The continuity relations are illustrated in Figure 8. The conditional probabilities satisfy a normalization 
property, which due to their piecewise form can be written as 

P+2,+ +P0,+ = 1 

P-2 - +P0- +P+2- = 1 

go (a;) + qi{x) + q2{x) = 1 

P-2,-+n{x)+r2{x) = l (17) 



As a result, one has to consider only 4 constants and 3 functions. The latter are solutions of coupled inhomogeneous 

nonlinear integral equations (more details are provided in Appendix B). We have solved the set of equations (B1-B2) 
by an iterative method, starting with p-2,- = 1, qoix) = 1, P-2,+ = 1, and all other functions equal to zero as initial 
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H 



FIG. 9: (Color on line) Average conditional probabilities as a function of the applied field for 2 = 3 (strong-disorder regime, 
(7 = 1.1): (a) p+2,+ and p_2,-, (b) po,+ and po,-, (c) P-2.+ and p+2,- (note the changes of scale). All probabilities are 
continuous functions of H\ p+2,+ , P-2,-, P-2,+ and p+2,- vary monotonically whereas po,+ and po,- go through a ma^ximum. 
The lines represent the solution of the exact self-consistent equations and the symbols the simulation results {Naitea ~ 10^, 

^realizationa 2.10 ). 



conditions. The "new" values of the PsgixY^ are then obtained by computing the integrals in the right-hand side of 
Eqs. (B1-B2) with the "old" values, using Simpson's rule. Convergence is assumed when maxjp^^™ — < 10~^. 
It is important to notice that, for a given external field, these equations may have several solutions (as discussed just 
below). When this happens, we find that only two solutions are stable, in the sense that they can be reached by the 
iteration algorithm. In addition to the solution obtained with the above mentioned initial conditions, the other one 
can be found by starting with p+2,- = Ij ^2(2^) = 1, <l2{x) = 1, P+2,+ = 1; and the other functions equal to zero. 
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FIG. 10: (Color on line) Average conditional probabilities as a function of the applied field for 2 = 3 (weak-disorder regime, 
(7 = 0.9): (a) p+2,+ and p-2,-, (b) po,+ and po,-, (c) P-2,+ and P-2,- (note the changes of scale). The probabilities show a 
singular behavior, with jump discontinuities, at nontrivial values of the applied field. The lines represent the solution of the 
exact self-consistent equations (both the physical and unphysical branches are shown) and the symbols represent the simulation 



results (Na-. 



= 10\ A^, 



realizationt 



= 2.10 ). A small finite-size effect is still visible in the simulation data in the jump regions. 



^From the solution of the coupled set of equations, it is straightforward to compute the average probabilities pg g{x). 
We find that the behavior drastically changes at a critical value of the random-field strength a. At large disorder, 
the set of equations has a unique solution and the average conditional probabilities Pg g are continuous functions of 
the applied field, as illustrated in Figure 9. At weak disorder, several branches of solutions are present. As explained 
above, two of them are "stable" . One corresponds to the physical process in which the magnetization always increases 
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with the apphed field H: the average conditional probabilities show a jump discontinuity as one increases H as 
illustrated in Fig. 10. The "unphysical" solution lies in a small range of fields. This solution is expected to be 
connected to the lower branch of the physical solution, but this intermediate portion is unreachable by the iterati! 
on algorithm (hence, is not shown in Fig. 10). We recall that this only concerns one side of the hysteresis loop, the 
other one being obtained by symmetry. 

This behavior is the signature of a disorder-induced out-of-equilibrium transition. Indeed, the magnetization which 
is given by Eqs. (7) and (8) shows a discontinuous jump as the applied field is increased only when the conditional 
probabilities are themselves discontinuous (for a value of x different form — J, J or 3J). The critical value of the 
disorder separating the two regimes, with and without a discontinuity, is equal to Uc ~ 1.006, as illustrated for 
P+2,+ {H) in Figure 11. 
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FIG. 11: (Color on line) Variation of the average probability p+2,+ (^^) with disorder strength a for z = 3 (physical branch). 
The critical disorder strength is found to be CTc — 1.006. The corresponding critical value of the applied field is He — 0.466. 



VI. DISCUSSION OF THE RESULTS AND COMPARISON WITH SIMULATION 

Before discussing the above results in more detail, we present the simulation study that we have carried out in 
parallel. Following the now standard procedure 9', .18], we have mimicked a Bethe lattice by random graphs of fixed 
connectivity. These graphs are easily handled in simulations because, even with a finite size, they have no boundary 
effects, contrary to the Bethe-lattice/Cayley-tree approach. However, the price to pay is that the random graphs have 
some loops. The (expected) equivalence between the two approaches stems from the fact that the typical loops found 
in random graphs have a length that scales like the logarithm of the total number of sites. 

The simulations are performed on systems of size NsUes — 10"* — 10^. For a Bethe lattice, such sizes are large 
enough to avoid significant finite-size effects and the results may be considered to be very close to the thermodynamic 
limit '17*1. In order to compute the magnetization as a function of the external field H, the latter is changed by a fixed 
amount AH — ±0.01. For a given realization {hi} of the random fields and a given value of H, we relax the system 
until all pairs are 2-spin-fiip stable, and we then compute the magnetization m ~ X^i^i ^i/^ ■ Similarly, to compute 
the conditional probabilities, we relax the system under the appropriate constraints. Finally, the results are averaged 
over Nreaiizations — 10^ — 10^ different realizations of random graphs and random fields. 

A comparison of the exact and simulation results is shown in Figures 6 (average conditional probabilities) and 7 
(magnetization) for the one-dimensional chain (z = 2) as well as in Figures 9 and 10 (average conditional probabilities) 
for the case z = 3. The agreement between simulation and theory is excellent. In the case z = 3 we actually exploit 
the agreement found for the average conditional probabilities to skip the analytical calculation of the magnetization 
(a feasible but very tedious task), and discuss instead the simulation results. 



The infiuence of the dynamics on the magnetization curve m{H) is illustrated in Figures 12 and 13. In Figure 12, 
we display the exact results obtained with the 1-spin-fiip and 2-spin-flip dynamics for the one-dimensional chain. The 
width of the hysteresis loop decreases as one goes from 1-spin fiips to 2-spin flips and it vanishes for the equilibrium 
curve (this latter is obtained by solving the exact recursion equations given in Ref.T^). A similar effect is observed 
for z = 3 : see Figure 13 where we have also added the curve obtained with a 3-spin-flip dynamics (whose precise 
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FIG. 12: (Color on line) Hysteresis loop, m{H), for the one-dimensional chain with the 1- and 2-spin- flip dynamics. The lines 
correspond to the exact results and the symbols to the simulation results (Nsues ~ 10^, NreaUzations = 10^). The dashed line 
is the equilibrium magnetization curve. The disorder strength is cr = 1.0. 

definition will be given below). This effect is quite understable from a physical point of view. Introducing more 
cooperativity in the dynamics allows the system to avoid being trapped in high-energy minima that are only stable 
with respect to single spin-flips. As a result, a better "equilibration" takes place and hysteresis is reduced. 
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FIG. 13: (Color on line) Hysteresis loop, m{H), for the z = 3 Bethe lattice with the 1- , 2-, and 3-spin-flip dynamics, as well 
as equilibrium magnetization curve (dashed line): simulation results for a — 1.0 {NsUes = 10^, NreaUzations = 10^). Note the 
jumps for the 2-, 3-spin-flip and equilibrium cases and the continuous behavior in the 1-spin-flip case. 



An interesting outcome of the present study is that the peculiar behavior of the hysteresis on a Bethe lattice with 
connectivity z = 3 found with the 1-spin-fiip dynamics disappears as one considers more cooperative dynamical rules. 
As discussed in the previous section, an out-of-equilibrium transition is indeed found with the 2-spin-flip dynamics. 
One can see from Figure 13 that not only does the hysteresis loop shrink as the dynamics is more cooperative, but 
it also becomes steeper; for the value of the disorder strength considered, it shows a discontinuity for the 2- and 
3-spin-flip dynamics (and for the equilibrium curve), whereas it stays continuous for the 1-spin-flip dynamics. 

The critical values of the disorder separating continuous from discontinuous behavior satisfy the following inequal- 
ities, 

^ - < af ^ ~ 1.006 < a^^ ~ 1.05, (18) 

where the superscripts \SF ^ 2SF, and EQ refer to the 1-, 2-spin-flip dynamics and to the equilibrium situation {<jf^ 
is estimated numerically by solving the recursion equations of Ref.|l3 for different values of the disorder strength). A 
similar set of inequalities was found in the computer simulation study of the RFIM on a 3-dimensional cubic lattice, 
with however a nonzero value of cr^^^ 7]. This suggests that increasing the cooperativity in the dynamics of the driven 
RFIM at r = 0, on top of reducing hysteresis effects, favors large-scale collective behavior of the spins. Accordingly, 
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the critical value of the disorder above which long-range correlations disappear increases with cooperativity. The 
absence of transition observed with the 1-spin-flip dynamics for z — 3 does not appear of physical significance. 

It is interesting at this point to discuss in more detail the way to increase the cooperativity in the spin dynamics. 
The idea is to introduce a /c-spin-flip dynamics 0, ^3 > i'^ which, as already seen for the 1- and 2-spin-flip dynamics, 
fc-stable states are visited. The main idea is that the system jumps from one fc-stable state to the nearest one whenever 
the former becomes unstable to the simultaneous flip of any subset of k or less spins. For the one-dimensional chain, 
implementing the rule is rather straighforward: one has to check the stability of continous strings of k' neighboring 
spins with k' < k. An exact description of the hysteresis behavior can again be derived. The relevant conditional 
probabilities now involve the configuration of a string formed by a spin in a given local field and the k — 1 spins 
Si, Sz-i at the lower levels, given that the parent spin is down and given the states of the k — 1 spins that are the 
descendan! ts of spin S^-i- The recursion relation that generalizes Eq. (9) requires transition rates for the change 
of configuration of a string of k spins, knowing the states of the parents and of the descendants. Writting down the 
equations is formally easy, the difficulty coming with the derivation of the explicit expressions for the transition rates. 
We therefore do not pursue here in this direction. Nonetheless, we point out that in the limit where k oo (k goes 
as the size of the system, here in one dimension) , one should likely recover the recursion expressions for the ground 
state of the RFIM 19]. Indeed, the /c-stable states explored by the system converge now to the ground state|l3. Il^ 
(expected to be unique for a Gaussian distribution of the random fields). 

In higher dimension or for Bethe lattices of higher connectivity, one has to worry about the geometry of the subset 
of k' spins. In the case of the 3-spin-flip dynamics on a Bethe lattice of connectivity z = 3, a simple implementation 
is to restrict the "geometry" of the triplets of spins considered (all single spins and all pairs of spins being nonetheless 
checked so that the 3-stable states so obtained are indeed 1- and 2-stable also): we have considered triplets formed 
by a spin and its two direct descendants, without checking triplets formed by a spin, a "child" in a given branch and 
a "grand-child" at one more level away (see Figure 2). Our simulations have been performed according to this rule. 



VII. CONCLUSION 



In this work, we have obtained the exact solution for the hysteresis behavior of the T = RFIM on a Bethe 
lattice when driven by an adiabatically changing applied field and evolved under a 2-spin-flip dynamics. This solution 
generalizes the results already derived for the 1-spin-flip dynamics. Although we have focused on the main hysteresis 
loop, it is most likely that further exact results are obtainable for other quantities such as minor hysteresis loops and 
the avalanche size distribution, albeit with considerable technical difficulty. 

Our exact results confirm the robustness of the whole picture concerning the nonequilibrium response of the RFIM 
under a change of spin dynamics robustness that was first shown by a simulation of the RFIM on a cubic lattice 0]. 
We find that for lattices with connectivity z > 3 the driven dynamics is characterized by an out-of-equilibrium critical 
point separating a strong-disorder regime with a smooth hysteresis curve and microscopic avalanches from a weak- 
disorder regime with a discontinuity in the hysteresis curve associated to a macroscopic avalanche. The width of the 
main hysteresis loop and the critical value of the disorder depend on the degree of cooperativity of the dynamics, but 
the overall scenario is unchanged. (As one expects the critical behavior on Bethe lattices to be mean-field-like, the 
critical point in the 2-spin-fiip dynamics trivially falls in the same universality class as that of the 1-spin-flip dynamics 
and that of the equ! ilibrium transition.) 

The robustness of the main characteristics associated with hysteresis behavior with respect to the spin dynamics 
suggests a nontrivial organization of the energy landscape. It was shown for the 1-dimensional chain and conjectured 
for the general case that the main hysteresis loop of the RFIM under 1-spin-flip dynamics is the outer envelope of 
all the 1-stable states in the magnetization versus applied-field diagramj23. A similar calculation could probably 
be undertaken for the 2-spin-flip dynamics, but in the absence of such a result, we merely speculate that the main 
hysteresis loop under /c-spin-flip dynamics is the envelope of all the fc-stable states. On physical ground, and as 
already shown here for the 2- and 3-spin-flip dynamics, one expects the width of the loop to shrink as the dynamics 
becomes more cooperative. One can then picture the organization of the metastable stable states as formed by a 
succession of envelopes, the kth envelope being included! in the {k — l)th and so on, until hysteresis vanishes and 
the envelope collapses to a curve that coincides with the equilibrium (ground state) behavior [22. As discussed in 
Ref. |2Cl| ■ a jump discontinuity associated with a macroscopic avalanche occurs when there is a reentrant part in the 
corresponding envelope: a jump in the fc-spin-flip dynamics would then imply a jump in more cooperative dynamics, 
but not necessarily in less cooperative ones. The apparently surprising result concerning the Bethe lattice with z = 3 
can be interpreted along theses lines: there is a transition for the 2-spin-flip dynamics, and for the 3-spin-flip dynamics 
and the equilibrium as well, but not for the 1-spin-flip dynamics. More work would be needed to confirm the general 
picture. 
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APPENDIX A: THE ONE-DIMENSIONAL CHAIN (z = 2) 

^Prom Eqs. (11), one can see that the conditional probabilities are constant in the intervals ] — oo, 0] and [2 J, +oo[, 

Pgg(x) = Pi;c^., whereas p^^{x) and p [x] have an explicit x-depcndence in the interval [0, 2J]. In this interval, the 

other probabilities are zero and the normalization property implies that pj^j^{x) = 1 — q{x) and p {x) = q{x), where 

q{x) is a function to be determined. Introducing the piecewise form of the conditional probabilities in the right-hand 
sides of Eqs. (11) leads to the following self- consistent equations: 

P++ =Pl-PO+ Po{p++ + P- + ) + (P2 - Pl)P++ 

P-+ = (1 - P2)P+- + (1 - P\)P— 

/■2J 

P+- = Po{p++ +P-+) + Sx {1- q{x)) 

Jo 

P— = {l-Pi){p+- +p — )+/ Sx q{x) 

Jo 

f—x p2J—x 

q{x) = {l-pi)p — +P+- Sxi+ 5xiq{xi) , (Al) 

Joo Jo 

where the p„'s are defined below Eq. (14). The normalization property in the intervals ] — oo, 0] and [2 J, +oo[ takes 
the form 

P+++P-+ = 1 

P+-+P— = l, (A2) 

which leads to further simplification of Eqs. (Al) that are reduced to 3 coupled linear equations, including one 
inhomogeneous integral equation for q{x). 

By using the definition of the average probabilities, Eq. (12), it is easy to rewrite the above equations in terms of 
those average probabilities and of q{x), and finally to derive Eqs. (13,14). 

Finding the expression for the magnetization, hence for the average probability that the central spin is up, requires 
the knowledge of the transition rates W^§_^ §^^^g,^ g,_^ _-j^^^^g^ g^ _^_-^^-^{xi,X2, x)) schematically pictured in Figure 4 (with 

only two branches merging at the central spin, since z = 2). The explicit expressions for the TVs are lengthy and 
not worth reproducing here, but can be derived in a straightforward manner by studying the stability of the triplet 
formed by the central spin and its neightbors under the 2-spin-flip dynamical rule. Inserting the expressions into Eq. 
(8) (with z = 2) and using Eqs. (12-14) lead to the following result: 

P{+l,\H)=p2{p++ +p-+f + 2pi{p++ +p-+)ip+- +p—) +Po{p+- +P—f 

+ 2-^[Fn(p++ +P-+) + F,o{p+- +P—)] - (/±^)'i^ii 
1 - Pi 1 - Pi 

+ 2(p++ -|-p_+) / 6x q{x) / 6xi +2{p+- + p — ) / Sx q{x) / 6xi 

Jo J -X Jo J2J-X 

^ I-2J I-2J pO j-2J px p2J 

- 2-^±^ / 5x q{x) / 5xi / 6x2-2 Sx q{x) / Sxi / 6x2 (A3) 

l~PlJo J2J-X J-xi Jo Jo J2J-X1 

where on top of the pn{Hys associated with 1-spin flips, we have introduced additional integrals over the (gaussian) 
distribution of the random field: the functions 

/.2J <.2J 

Foo{H) = 6x 6x1 
Jo J2J-X 

r2J pO 



p^J /■U 

F,o{H) = Foi{H) = / Sx / Sx, 

Jo J -X 





/O M 
Sx / Sxi (A4) 

-27 J-(2J+x) 

represent the probabilities that a pair of neighbors flips up in an irreducibly cooperative way, given that the neighbors 
of the pair are both down, one up and one down, and both up, respectively, whereas 

p2J pO 

Fn{H) = Sx{ Sxif (A5) 

Jo J-x 
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is related to the probability to flip a spin up by an irreducibly cooperative flip of any of the two pairs it belongs to, 
knowing that the outer neighbors of the pairs are up. Again, the dependence of all these quantities on the applied 
field H comes through the measure 6x (see below Eq. (6)). 

Notice that the first line of Eq. (A3) is exactly the result obtained with single-spin flips only 8] . One can check that 
by setting q{x) = (1 — ?'i)/(l — Pi + Po) over the whole interval [0, 2 J] and using Eqs. (13), the other contributions 
cancel altogether. 

The magnetization is finally obtained through Eq. (7). 



APPENDIX B: SELF-CONSISTENT EQUATIONS FOR THE CASE 2 = 3 



The transition rates M^|g|(5' 5/ „i)^(5j^ 5)|_i|(a;i, a;2, x) corresponding to the case 2: = 3 in Figure 3 are ob- 
tained by combining two stability diagrams similar to those of Figure 1 in a 3-dimensional graph with axes xi,X2 
and X. Deriving the expressions is a tedious task, but with no conceptual difficulty. The expressions lead to the 
piecewise form of the conditional probabilities given in Eqs. (16) with 4 different intervals of x to be considered: 
] — 00, —J], [—J, J], [J, 3J] and [3J, -|-oo[. The probabilities are independent of x in the intervals ] — 00,— J] and 
[3J, +oo[ with values pgg (see Eqs. (16)) given by the self-consistent equations 



P+2,+ 



Po. 



Po, 



P~2. 



p-JiS-l) POO 

V/ dxpg_{x)+y2 pg \ 

5 i-J(S+l) 5 i-J(S-l) 



x) 



= 2 



P-2.+ = 



P+2, 



E 
E 

s 

E 

s 

E 
E 



-J(S+1) 



7(S-1) 
7(S-1) 



Sx Ps-{x) 

Sx Pg _{X) 

Sx Ps,+ {x) 
Sx Ps.-(x) 



E 



-J(S-l) 



-J{S+1) 



Sx Ps,-W + E 



E 



J(S-l) 



Sx Ps.+{x) 



-J{S-1) 



Sx Pg _ (X) 



J(S-l) 



Sx Ps,+ {x) 



(Bl) 



Note that sums and integrals do not commute in the above equations. In the intermediate intervals [—J, J] and 
[J, 3J], the probabilites are functions of a; (except p-2,- m [—J, J]). The corresponding functions, introduced in Eqs. 
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(16), satisfy the following equations: 

»-,7(S-l) 



+ 



qi{x) = 2 



90 (a;) 



r2{x) 



E 

5 

Y. 

E 

§ 

E 

s 

E 



SXlPg _{Xl) 

J(S-2)-x 

X} 

5xlp§_^_{xl) 



J(S-l) 

J{S+1) 



SxiPg_{Xi) 



-J(S-l) 
J{S+1) 
-J(S-l) 

J(S+1) 
-J(S-l) 



5xip§_{xi) + 



-Jis-i) 



^E / ^xip§_{xi)-^ I ^ Sxipg_{xi) 

s 

Sxipg {xi) + Y, _ ^xiPs 

-J(S-2)-x ^J-J(S-l) 



J{S-2)-x 

SxiPg_^{xi) 



J(S-l) 



{xi) 



-J{S-2)-x 



-J(S-l) 
JS-x 
-JS-x 



Sxipg _{xi) 



^XiPg_{Xi) 



J(S-l) 



<5a;iPs,+(a;i) 



J(S-l) 



6xipg_^_{xi) 



5xipg_{xi) 



E 



J(S-l) 



5xiPs,^{xi) 



(B2) 



Note that the x-dependence appears through the bounds of the integrals. In Eqs. (Bl) and (B2) the conditional 
probabilities present in the right-hand sides should be replaced by their expressions given by Eqs. (16), which makes 
the whole set of equations self-consistent. Due to the normalization conditions there arc only 7 independent coupled 
equations. All the equations are nonlinear, and those for the 3 independent functions are inhomogeneous integrals 
equations as well. The procedure to solve the set of equations is discussed in the main text. 

Finally, we point out that, as in the z = 2 case, the result of the 1-spin-flip dynamics can be recovered by assuming 
that the functions in the intermediate range [—J, J] and [J, 3 J] are constant. The proof is a little more complicated 
than for 2; = 2 and is not worth reproducing here. 
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